Theory and applications of the stress density 
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Drawing on the theory of quantum mechanical stress, we introduce the stress density in density 
functional theory, and give specific prescriptions for its practical and efficient implementation in the 
plane wave ultrasoft pseudopotential method within the local-density approximation. In analogy 
with the Chetty-Martin energy density, the stress density provides a spatial resolution of the contri- 
butions to the integrated macroscopic stress tensor. While this resolution is inherently non-unique 
(gauge-dependent), there exist gauge- independent ways of using it in practice. Here we adopt the 
following ones: a) calculating integrated macroscopic stresses over appropriately defined parts of a 
system; b) analyzing macroscopic averages of the stress density; c) analyzing changes in the stress 
density in response to external perturbation. The abilities of the stress density are demonstrated 
for a set of representative test cases from surface and interface physics: in perspective, the stress 
density emerges as vastly more powerful and predictive than the integrated macroscopic stress. 
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I. INTRODUCTION 

Energy and stress are quantities of obvious fundamen- 
tal importance to the physics of solid state systems. In 
the last fifteen to twenty years, it has been possible to cal- 
culate them with progressively increasing accuracy from 
first-principles using density- functional theory (DFT)£I 
The energy and stress as output from a typical calcula- 
tion are a scalar and, respectively, a tensor whose val- 
ues pertain to an entire finite system or to the unit cell 
of a periodic system. The relevance of the information 
they carry is concealed in their parametric dependence 
on external constraints, such as atomic geometries. It 
is not difficult, on the other hand, to picture situations 
in which access to local information related to energy 
or stress would be highly desirable. Clearly, such access 
would require a procedure allowing to sort out contribu- 
tions to energy or stress due to specific regions inside the 
simulation cell. 

By way of example, consider the calculation of the en- 
ergy or stress change due to surface formation by means 
of the common repeated-slab supercell technique. The 
energy (or stress) of the slab contain several informations 
coupled to one another, which one would like to access 
separately. A simple case is that of a symmetric slab: to 
extract the properly surface-related quantities, the bulk 
contributions must be subtracted out, and this requires 
an independent bulk calculation. A qualitatively more 
complex case is that of an inherently asymmetric slab, 
such as frequently found in polar semiconductor surface 
simulations: in that case, two inequivalent surfaces are 
present, and their individual surface energies cannot be 
calculated independently. [It is occasionally possible to 
symmetrize the simulation cells in such a way as to obtain 
the surface energy of one of the surfaces, but even then 
only against a major computing effort.] Clearly, more 



basic issues may be addressed by a local energy or stress 
analysis: if a surface under stress reconstructs, and stress 
is the alleged driving force, it would be enlightening to 
know which parts of the systems tend to reduce rather 
than increase their stress during the transformation. 

Chetty and Martin,El looking for conceptual as well as 
practical solution to these and similar issues, introduced 
the energy density in the DFT framework. The energy 
density is a position-dependent density-functional scalar 
field whose integral upon the unit cell equals by construc- 
tion the total energy. In general, this quantity is gauge 
dependent, i.e. not uniquely denned. Indeed, due to the 
long-range Coulomb interactions, the correspondence be- 
tween a region of space and the energy stored therein is 
arbitrary, unless the value of the functional is fixed by 
symmetry at the boundaries of the region. 

In an isolated system, the energy density is zero at the 
boundaries and outside, and the integral over the region 
enclosed by the boundaries is unique, and equals by con- 
struction the total energy. A more interesting situation 
where these requirements are satisfied is a surface (or in- 
terface) slab, infinite and periodic in the surface plane, 
and finite in the orthogonal direction. If its thickness is 
sufficient that bulk-like behavior is recovered in its inte- 
rior, the energy associated with regions bounded by bulk 
layers is uniquely fixed and equal to the bulk .energy (also 
relevant to the simple operative definitions^ of surface 
energy). The vacuum contribution is zero by definition. 
Thus, the energy density in the surface regions obeys the 
appropriate boundary conditions, hence its integral upon 
that region is uniquely determined and gauge invariant. 
Similar arguments apply to interface slabs, that couple 
two different bulk regions and one (or quite possibly two, 
independent) interface region(s). 

The concept of quantum-mechanical microscopic stress, 
was introduced in the early days of quantum field theorycl 
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and revisited and deeply investigated more recently.ETEl 
In this paper, building upon the modern microscopic 
quantum-mechanical theory of stressaH and drawing on 
the Chetty-Martin theory of the energy densityp we con- 
struct a position-dependent stress density. We start by 
defining the stress density as a rank-2 cartesian tensor 
field T a p(r) whose integral upon the cell equals by con- 
struction the macroscopic stress: 



T. 



a/3 



drT a/3 (r) 



(1) 



The formulation we present below is adapted to the prac- 
tical plane-wave pseudopotential computational frame- 
work. The alternative definition of a stress field related 
to the force field Fp(r) through 



f>(r) = -V a S af) (r), 



(2) 



cannot be used in this framework, because it fixes S 
only up to a gauge (the curl of an arbitrary function). 
While the gauge cbpice is unrestricted in principle, only 
the Maxwell gaugeQ is practically tractable for the stress 
field: that gauge satisfies Eq. (0) only for systems 
with purely coulombic interactions, conflicting with the 
fact that pseudopotentials produce non-coulombic inter- 
actions. 

Numerous applications can be envisaged for the stress 
density. Those we focus on here concern surface and 
interface physics, as indeed structural instabilities and 
specific electronic effects at surfaces are fairly commonly, 
associated with (among others) the surface excess stress.El 
Three point are worth some attention. First, a spatially- 
resolved stress function can be used to disentangle the 
macroscopic, integrated stress of inequivalent surfaces (or 
interfaces) that may be coupled into a single simulation 
cell, as is frequently the case for semiconductor polar 
orientations. Second, residual bulk contributions to the 
integrated stress can be computed from the stress density 
in the same calculation in a technically consistent fashion. 
Third, but definitely not least, the stress density can be 
applied as a "stress microscope" to investigate the role of 
stress in specific surface or interface processes; this kind 
of analysis, not relying merely on integrated macroscopic 
stress, can change qualitatively our understanding of. the. 
relation of surface stress with (e.g.) reconstructionErEJ 
and adsorption. 

The plan of the work is as follows: in Sec. |FJ we for- 
mulate the stress density for density-functional theory in 
the local density approximation (LDA), for the specific 
instance of the plane- wave pseudopotential framework. 
In Sec. Ill we test the functionalities and accuracy of 



the stress density: specifically, in Sec. Ill A we apply it 
to the Al (111) surface, in Sec. Ill B to the t he po lar 
(100) and (111) surfaces of GaAs, and in Sec. Ill C to 
the ( 001) and (110) Ge/GaAs interfaces. Finally in Sec. 
Ill D| we discuss an application of the stress density as 
a microscopic analysis tool in the case of metal surface 
reconstructions. Atomic units are used throughout. 



II. STRESS DENSITY IN DFT-LDA 
A. Detailed formulation 

Nielsen and MartirH carried out the quantum formula- 
tion of macroscopic stress for a general system of interact- 
ing particles. They also presented!] an expression suitable 
for the DFT-LDA framework, and the explicit form of the 
stress tensor for a plane^wave basis and pseudopotentials 
(a generalization existsEJ to gradient-corrected exchange- 
correlation functionals) . An analogous focmulation for an 
LCAO basis was provided by Feibelman.EJ 

The Nielsen-Martin expression for the macroscopic 
stress is comprised of kinetic, exchange-correlation, 
electron-electron interaction, and electron-ion interaction 
term. For each of these contributions, we derive a space- 
resolved expression whose integral over the simulation 
cell equals the corresponding macroscopic quantity. 

The macroscopic kinetic stress can be expressed in two 
equivalent ways, which reflect the gauge dependence of 
the kinetic stress density. The symmetric form is 

= e 2 J2 /xkWk / drV a 4(r)V^ k (r), (3) 

and the antisymmetric form is 

T «P = ~ e2 E ^ kCJk / dv €k(r)V Q V^,k(r), (4) 



//k 



where and lu^ are occupation numbers and k-point 
weights, respectively. The symmetric expression in Eq. 
@ is preferable computation-wise as it does not involve 
second derivatives of the wayefunctions (for other reasons 
to prefer this form, see Ref .□) . Thus we define the kinetic 
stress density as 

7"i1 n (r) = %r /^k[V a Ck(r)V^,k(r) + c.c.]. (5) 



The calculation of 7^ n (r) is computationally inexpen- 
sive. The wavefunction derivatives are promptly calcu- 
lated in Fourier space and then carried back to real space 
by fast-Fourier transform (FFT). 

The exchange-correlation contribution to the macro- 
scopic stress is 



T% = -&«(> / dvp{v)[e xc {v)-V xc {v)l 



(6) 



where e xc and V xc are the exchange-correlation energy 
density and the exchange-correlation potential, respec- 
tively, and p is the electronic charge density. The 
corresponding expression for the stress density follows 
straightforwardly: 



2 



T^(r) = -p(r)[V xc (r) - e xc (r)]S aP . 
The macroscopic Hartree stress is given by 

,,2 



(7) 



e 

rpG—e __ 

J ^ - 2 



/|** V(r>(r0 A_l_ 



= i d rp(r)p(r } — ' (8) 

Due to the long-rangedness of coulombic interactions, 
there is no unique way to fix the stress density as a func- 
tional of the charge density. However, an optimal choice 
is the Maxwell gaugejdu which is physically transparent 
and computationally easy to implement in our specific 
framework. The Maxwell stress density is 

T «» = -^[E a {v)E p {v) - \5 a0 E\v% (9) 



where E n is the electric field 



E a (r) = -V Q / dr 



J Pi*') 



(10) 



To cure the coulombic divergences in 7^9 (r), the total 
charge density 

should be used in Eq. (|l0|), rather than the electronic 
charge density p. The ionic charge is represented by a 
sum of ion-centered Gaussians as 



.3/2 £3 



(11) 



where Rj and Zj are ionic positions and charges, re- 
spectively, and R c the gaussian radius. In this way 7^ 
contains the electron-electron interaction term, 



plus an unphysical term due to the fictitious interaction 
between electronic and Gaussian-ion charge density, to 
be taken care of below (see Eq. (|l7|)). This procedure 
ensures charge neutrality within the volume bounded by 
bulk layers, besides that of the whole simulation cell. The 
crucial requirement that the stress density integral over 
a bulk layer inside, e.g., a surface slab to be equal to 
the bulk stress, is then satisfied. The electric field can 
be calculated in G-space and then carried to real space 
where T^(r) is easily evaluated. 

A further contribution to the ion-ion interaction arises 
from the strain derivative of the ionic charge density. In 
real space, this Ewald-like term reads 



^-Ewald /„\ 

1 a/3 \ V ) — 



r>2 



'(r)4°"(r) 



(14) 



For a purely coulombic system, the previous formu- 
las are all that is needed to construct the stress density. 
However, a practical implementation for the plane-wave 
pseudopotential method calls for several extensions, due 
mostly to the non-coulombic contributions of pseudopo- 
tential. The local pseudopotential stress contribution is 



rpC — i 

J q/3 — 



"drp(r)^ — V}° c 



'(r-R,), (15) 



where Vj° c is the local part of pseudopotential sited on 
atom j. From d\r\/de a f3 — r Q r^/|r|, it follows that 

t ( r - Rj) a (r - Rj)f3 

i J 

(16) 



where V' l ° c (x) = dVj° c (x) /dx. We define the corre- 
sponding stress density as 

d(V} oc (r) + V} on (r)) 
(r) = p(r) £ J [ > 3 (17) 



1 [Et%v)Ey{v)- l -8 aP E^ 2 {v)l (12) 



and the ion-ion interaction 



^{E™(r)E™(v) - h af3 E™ 2 (r)], (13) 



where the potential Vj° n generated by the the Gaussian 
ionic charges is added to compensate the same contri- 
bution present in Eq. (^]), and to cancel the coulombic 
divergences in T c ~\ The strain derivative of local pseu- 
dopotential in G-space is (for brevity we use the same 
symbol for a quantity in G- and r-space) 



BV loc 

^-(G)=5>(G') 



dV' oc {G') Aire 



dG' 



^±-pf»(G>)) 2 G' a G' + (v}° c (G>) + Iprpf 1 ^)) <W 



(18) 



for all G^O, whereas the G = contribution is 



gyjoc 

de a0 



(0) 



V/3 



(19) 



Eq. (|l|) is derived using the fact that a positive strain 
e in r-space corresponds to linear order to G' = (1 — e) G 
in G-space, thus 
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<9V; loc (G) dV} oc (G) 

The G-space derivative of Vj oc is evaluated numerically. 

The non local pseudopotential term only acts in the 
core regions, where spatial resolution is largely arbitrary. 
Therefore, we choose it to be a superposition of ionic-site 
centered contributions, 

^( r )=E P a^( r - T ^ ( 21 ) 

i 

where stress 

gyNL 

Kb = Y,fukVk{il>vk\-z*—\rl>uk) (22) 



is the j -ionic site contribution. The most general form 
of a fully non local pseudopotential can be written as a 
non-diagonal projector 

Vf L = J2D J nm \(3i)(f3L\, (23) 

nm 



where n, m are sets of atomic quantum numbers. In 
the case of ordinary Kleinman-Bylander pseudopoten- 
tials, the matrix D 3 nm is diagonal and its elements are 
constant, thus the corresponding stress is given by the 
strain derivatives of the projector functions (3 3 n . In plane 
waves, we have in general 



r aB — r aB 



k 



D 3 t 

nm 



-i(G-G')-R 



' Ck(G) 



d[P 3 n (k + G)/3 3 m *(k + G')} 
de a B 



^k(G'), 



(24) 



where 

»k + G) 



Be, 



aB 



The term P; 



3/%(k + G) 
d{k + G) a 



{k + G)g. 



(25) 



is only non-zero if ultrasoft 



pseudopotentialso are used. In that case, D 3 nri 
diagonal and contains a charge-dependent term, 



is non 



Dim = D°nm + J *T <&Jt - RjW^T) 



+ F Hxc (r)] 



(26) 



where ]/ Hxc is the screening (Hartree plus exchange- 
correlation) potential. The strain derivative of |V loc (r) + 
I/ Hxc (r)] was already taken into account in Eq. (|l7|), so 
the strain derivative of D 3 nm only contributes the term 



r aB 



with 



V B 3 



-iG-R, 



(V 



loc 



V Hxc )*(G 



dQ nm (G) 



G 



dc 



aB 



B 3 



/«* Wk(^k|/%)(/%JVVk>. 



(27) 



The formulation of the stress density is now completed. 
The final expression reads 



T a0 (r) = T$Xr) + (r) + 7?S(r) 



(28) 



1 aB 



L " iad (r) + r Q y(r) + r^(r), 



macroscopic stress tensor. In order to accurately evaluate 
the stress density and the stress itself, a well-converged 
selfconsistent charge density is generally needed, due to 
the non-variational character of the stress, and its sensi- 
tivity to the details of the charge density. 



B. Averages and stress density handling 



A practical point concerns how to handle and visual- 
ize the stress density (more discussion and examples are 
given in Sec. III). When dealing with surfaces and inter- 
faces, appropriate averages of the stress density can be 
taken, that contain all the relevant information. If z is 
the surface normal, the planar average is defined as 



the involved terms being given respectively by Eqs. (|^), 
(@), ®> (@> and ©■ Operatively, all terms (ex- 
cept exchange-correlationjare first calculated in G-space 
and then Fourier-transformed to r-space, with a compu- 
tational load comparable to that needed to calculate the 



(29) 



iG,-z 



A typical planar average (see Fig. [I]) is strongly oscil- 
lating, with negative peaks at the ionic positions (corre- 
sponding mostly to attractive ion-electron contributions) 
and positive peaks in the interstitial region (correspond- 
ing mostly to kinetic compressive stress). Deviations 
from bulk- like behavior, signaling surface- induced stress 
changes, can be directly inspected. 

A further convenient way of analyzing T(r) is its 
macrosopic averaged over a period of length d, and de- 
fined as 



T{z) 



1 



z+d/2 



dz' T(z') 



(30) 



z-d/2 
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This filtering operation eliminates oscillations of period 
d, and can be generalized!!!] to account for several super- 
imposed oscillating components. If d is, e.g., the inter- 
layer distance in a metal slab, then T(z) is constant in 
the bulk region of the slab, and its value equals the bulk 
stress; in the vacuum region T vanishes by construction. 
Any deviation from these constant values in the simula- 
tion cell quantify surface- related stress changes. In G- 
space the macroscopic average reads 



T(G Z 



2 no* 

d G z 



dG z 



(31) 



where T{G Z ) — T(0,0, G z ). As discussed in detail in 
Ref. ||, the macroscopic average is gauge-independent by 
construction. 



C. Gauge dependence 

A major concern about the applicability of the stress 
density (as in the energy density case) is its formal gauge 
dependence, showing up e.g. in the two distinct forms 
Eq. (||) and (||) of the kinetic stress contribution. Gauge 
dependence, however, does not affect the three classes of 
applications we are interested in here, which involve in- 
tegrated stresses, macroscopic averages of the stress den- 
sity, and changes in the planar averages thereof. 



1. 



By definition, the gauge indeterminacy does not af- 
fect observables such as the integrated stress, hence 
it does not prevent the usage of the stress density 
to extract, say, integrated surface stresse s in an im- 
proved or technically safer fashion (Sec. 



irj 



2. As discussed in detail in Ref. ^J, the macroscopic 
average of the stress density is gauge-independent 
by construction, since gauge functions are lattice- 
periodic, and hence are filtered out by macroscopic 
averaging. 

3. Changes in stress density in response to perturba- 
tion are gauge-independent. 

The proof of the last statement is as follows. Let T[p(r)] 
and T'[p(r)] be two functionals of the density whose in- 
tegral is the same: 



dvT[p{v)]= J dvT'[p(v)}. 



(32) 



[An interesting choice is the two functionals being the 
stress densities with the two kinetic gauges Eqs. (gj) and 
(0).] Under a generic perturbation on the Hamiltonian, 



H X = H + XV, 
T becomes to linear order 



(33) 



r = t[ Po ) 



d Pj P0 

As the unperturbed integrated values are equal, 
J drT[ Po ] = J dvT'[p }, 



(34) 



(35) 



substituting Eq. (|34|) into Eq.(|32|) we get 
= Jdr (T[p }-T'[p ]) = 



drSp x 



dT' \ _ fdT\ 



dp 



d pJo, 



(36) 



The density variations are arbitrary, so we conclude that 



9T' 

~dp~ 



(37) 



Therefore, two density-functionals obeying Eq. £32] have 
identical variations under external perturbation, to lin- 
ear order with respect to the density. We then rest as- 
sured that stress density changes are gauge-independent 
to linear order. This puts on fir m grou nd the kind of 
application exemplified in Section HID, which rely on 
the analysis of changes in stress density planar averages; 
typically, one will compare e.g. planar averages in differ- 
ent systems and instances instead of analyzing differences 
(which is justified, as just proven): the reason for this is 
simply the better pictorial rendering of the former ap- 
proach. 

As to stress densities in themselves (more precisely any 
non-macroscopically averaged ones), they are unavoid- 
ably non-unique. In practice, however, we found numer- 
ically very similar results using both the possible choices 
for the kinetic stress density; in addition, the spatial be- 
havior of e.g. planar averages agre es closely with physical 
intuition (as discussed in Sec. III). 



III. APPLICATIONS OF THE STRESS DENSITY 

A. Non-polar surfaces: tests on Al (111) 

To illustrate the basic characteristics of the stress den- 
sity we first discuss the (111)— 1 x 1 surface of Al in the 
ideal (i.e. unrelaxed) structure. Calculations employed 
a 9-layers slab, norm-conserving pseudopotentials, and 
well converged plane-wave and k-point sets. In our for- 
malism a negative stress is tensile, i.e. it indicates that 
reducing the interatomic distance is energetically favor- 
able; a positive stress is compressive - that is, the system 
prefers to be stretched. The stress densities considered 
in the following are planar or macroscopic averages of 
the in-plane components T xx (r) = T yy (r) (almost all the 
systems considered hereafter are isotropic in the surface 
plane): indeed, only the in-plane stress matters, since the 
orthogonal component is eliminated by relaxations. 
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FIG. 1. Stress density of the unrelaxed Al (111) surface. 
Solid lines are planar and macroscopic averages for the sur- 
face slab, dashed lines same quantities for the bulk in the 
(111) slab symmetry (see discussion in text). 

In Fig. | four curves are presented: two of them, the 
solid lines, are the planar and macroscopic averages of the 
stress density for the surface; the other two, the dashed 
lines, are the same quantities for the Al bulk in (111) 
symmetry (obtained filling up with atoms the vacuum 
region between repeated slabs). The comparison serves 
to better visualize changes due to the surface formation. 
Thanks to symmetry along 111], only half of the simu- 
lation cell is shown in Fig. [l. 

As mentioned earlier, the planar averages are charac- 
terized by oscillations having the period of the interlayer 
distance in the bulk region. Negative peaks correspond to 
the layer positions by construction, and mostly store the 
tensile contribution due to the electron-ion interaction, 
whereas the positive peaks located between layers are 
dominated by compressive kinetic and coulombic terms. 
The surface planar average drops to zero on entering the 
vacuum region, and it is informative to compare it with 
the bulk planar average near the surface. 

The macroscopic average of the surface is constant in 
the bulk region and zero in vacuum,?. e. it is not constant 
only in the region where surface effects arise. Thus, it is 
enables one to identify the boundaries of the surface re- 
gion. In the present case, for instance, the perturbation 
produced by surface formation extends roughly over one 
interlayer distance at the solid-vacuum interface. The 
macroscopic average for the bulk is everywhere constant, 
and it equals the surface macroscopic average in the bulk 
region. This common value of the two macroscopic av- 
erages in the bulk region, equaling by construction the 
bulk stress, is close but not equal to zero. 

This is a first important capability of the stress den- 
sity. It is a common feature of stress calculations that the 
stress be non-zero even though the lattice constant cor- 
responds to theoretical equilibrium (hence in principle to 
zero stress). Indeed, energy minimization is usually per- 
formed at fixed plane- wave cutoff energy (i.e. variable 



number of plane- waves) , whereas the stress.ia implicitly 
calculated at fixed number of plane waves.E^EH This bias, 
specific of the plane- wave methodology, can be eliminated 
using an absolutely converged basis set, but that is highly 
impractical especially when dealing with large systems. 
Thus, the integrated surface stress should be obtained 
substracting n-times the bulk contribution (if n is the 
number of atoms, or more generally of formula units, in 
the slab) from the total stress as 



1 af3 — 2 L a 



lab 





(38) 



where factor 1/2 accounts for the two slab surfaces. As 
for the surface energies,tJ this procedure is justified in 
general to the extent that the finite-size basis error is 
proportional to the number of atoms, with the additional 
proviso that the error per atom does not change signifi- 
cantly for small changes in volume per atomJU 

The stress density provides a straightforward evalua- 
tion of the spurious bulk stress, requiring no separate 
calculations, because the spurious bulk stress per atom 
is just the value of the macroscopic stress density aver- 
age in the bulk region. Calculating the spurious stress 
in this fashion not only saves computing time, but also 
gives more accurate results, since use of the bulk stress 
calculated in the bulk unit cell in Eq. ( |38| ) introduces 
numerical errors beyond the basis-set error, due e.g. to 
inequivalent k-points. Another source of error is surface 
stress anisotropy (such as it occurs in most reconstruc- 
tions): the spurious bulk stress is also expected to be 
anisotropic, a feature that cannot be but missed in the 
bulk unit cell. 

Coming now to the shape of the macroscopic stress 
density average, we see that the surface formation causes 
an excess of tensile stress just below the surface- vacuum 
interface, partially compensated by an overshot of com- 
pressive stress outside the interface. The resulting total 
surface stress is slightly tensile. The presence of tensile 
stress is a common feature of most metal surf aces, E3~ — 
and its origin is easily understandable for sp metalst£ 
in terms of electron spill-out into vacuum, resulting 
in a reduction of the kinetic compressive stress in the 
near-surface region (the interpretation is slightly more 
involvedta for transition metals). This tensile stress re- 
mains built-in at the surface, since surface atoms are con- 
strained to their positions by the substrate potential, and 
in-plane contraction would require the formation of (typ- 
ically costly) defects to be topologically possible. 

We close this Section with an example of layer-by-layer 
resolution (LBL) of the stress. In Table | we list the LBL 
of the total stress for Al (111), i.e. the integrals of the 
planar stress density average over one interplanar dis- 
tance. The 15 "geometrical" slab layers (9 are occupied 
by atoms) are labeled -7 to 7, layer being the bulk- 
like slab center, and layers -l-n and -n being equivalent 
by inversion symmetry along z. Layers -3 to +3 are 
bulk-like, having nearly equal LBL value. Their aver- 
age stress (-0.743 eV) equals to within 1 meV the bulk 



G 



stress calculated in the bulk fee. Inserting this value into 
Eq. we obtain a value for the surface stress (-0.58 
eV/atom) unaffected by numerical errors. The resulting 
surface stress is in good agreement with the results of a 
similar procedure reported by FeibclmanJij 



B. Polar surfaces: GaAs (100) and (111) 

As briefly mentioned earlier, orientations for which 
cleavage produces inequivalent surfaces are a natural 
playground for the stress (and energy) density. We con- 
sider here two prototypical cases, namely the (100) and 
(111) surfaces of (zincblende) GaAs, which have-been 
and still ace. the focus of intense experimentalEirEfl and 
theoreticaOtil work. These surfaces are polar, and 
present several (otherwise unpleasant) features making 
them suitable for a test of the stress density. 

For the (100) orientation, each surface (Ga or As ter- 
minated) can indeed be individually simulated in a sym- 
metric slab, so that the energy of each surface can be 
calculated via total energy differences; for testing pur- 
poses, however, we will also consider asymmetric slabs 
where inequivalent surfaces are coupled. For the (111) 
orientation, the surfaces are intrinsically different. Ei- 
ther are they geometrically identical but terminated by 
a different atom, or terminated with the same atomic 
species but with different geometries: the (say) Ga- 
terminated surfaces of a (lll)-oriented slab are neces- 
sarily one (111) and one (111), and no symmetrization is 
practicable inJ.bis.case. Using the energy density is then 
indispensablec3~E3 for a direct calculation of the abso- 
lute surface energy of each surface, whereas only aver- 
ages, or differences to some reference system are accessi- 
ble by plain total energy calculations. The procedure is 
not straightforward, since it requires counting the atoms 
that belong to a given surface; this can be troublesome 
for GaAs 1), and arguments based on cell symmetry 
adaptations or space partitioning^ have to be employed 
to complement the use of the energy density. Neverthe- 
less, the energy density is more computationally conve- 
nient to study this surface than the total energy is. 

Analogous considerations hold for the stress density, 
with the significant difference that there is no need to 
count surface atoms to evaluate the surface stress, and 
one need only making sure that the bulk structure is 
at theoretical equilibrium and the cutoff energy is large 
enough to render negligible the finite-size basis error. 

There are other difficulties related to the presence of 
unequivalent surfaces. The two workfunctions are not 
lined up, and fictitious electric fields are produced in the 
vacuum region. Furthermore, as-cleaved GaAs (100) and 
(111) are metallic, because partially filled dangling-bond 
surface states are present on Ga- as well as As-terminated 
surfaces. As a consequence, a sizable charge transfer 
between the surface bands occurs across the slab; this 
artificial charge flow can be avoided resorting to passi- 



vation with, typically, fractionally charged H atoms. In 
reality these surfaces reconstruct into complex atomic ar- 
rangements to restore neutrality and remove metallicity; 
however, we decided to consider here only the unrecon- 
structed and unpassivated surfaces, since they represent 
a technically intricate and critical case of supercell sur- 
face calculation if there ever was one, and are ideally 
suited as a testing ground for the stress density. We left 
the surfaces unrelaxed at the theoretical lattice constant. 

Ga-terminated (Ga-t) and As-terminated (As-t) GaAs 
(100) can be represented individually within symmetric 
slabs. In this way, essentially no charge transfer across 
the slab occurs, barring the small residual interaction be- 
tween surfaces which splits the (in principle degenerate) 
surface bands. On the other hand, the Ga-t and As-t 
surfaces can also be made to appear simultaneously in 
the same cell. This configuration is undesirable in gen- 
eral as it leads to sizable charge transfer, but it is useful 
as a test case; surface stresses obtained in the two sym- 
metric (charge transfer unaffected) slabs can be directly 
compared with those of the coupled-surfaces asymmetric 
slab. For these lxl surfaces we use 12 (plus 6 vac- 
uum) layers for the asymmetric slab, 11 (plus 7) layers 
for the symmetric one. Ultrasoft pseudopotentialsE-3 are 
employed. The macroscopic averages are obtained by av- 
eraging over a distance of two interlayer spacings. 

In Fig. ^| we show planar and macroscopic averages of 
the stress density in the asymmetric slab (solid lines), and 
the same quantities in the symmetric Ga-t slab (dashed 
lines) . The macroscopic average inside the slab show that 
the bulk stress is zero to numerical accuracy, so that 
the surface stress is given directly by stress density in- 
tegrals over one half of the slab. [Because the surfaces 
are non-stochiometric, Eq. ( |38| ) cannot be used to elim- 
inate the spurious bulk stress; we therefore had to use 
a sufficiently converged (60 Ryd cutoff) plane-wave ba- 
sis so as to achieve a negligible residual bulk stress.] In 
the bulk region the averages for the symmetric and asym- 
metric slabs are indistinguishable, indicating that in both 
structure the inner slab region is insensitive (stress-wise) 
to the surface. At the Ga-t surface the various aver- 
ages match almost perfectly, indicating that the slabs 
are large enough to entirely decouple (stress-wise) the 
surfaces (mind that in one case the other surface is iden- 
tical, in the other it is As-t). As expected, the calculated 
charge transfer is much larger for the asymmetric (0.02 
electrons per cell, or 5 mC/m 2 ) than for the symmetric 
slab (0.5 mC/m 2 ). The surface stress per surface cell ob- 
tained integrating the stress density upon the half slab 
containing the Ga-t surface is -1.114 eV/atom per surface 
cell for the asymmetric slab, and -1.148 eV/atom for the 
symmetric one. Thus the error due to charge transfer is 
about 3%, and well below 0.1 eV per surface cell, which is 
similar to that of half-slab integrals of the energy density 
in Ref. where an analogous test was performed. Notice 
that for these unrelaxed structures the surface effects are 
very much localized on the surface atoms, and the bulk 
features are substantially recovered just one layer below 
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the surface. Also, not surprisingly, the surface stress is 
tensile. 
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FIG. 2. Planar and macroscopic averages of the stress den- 
sity for GaAs (100). Solid lines refer to the (asymmetric) slab 
containing both Ga-t and As-t surfaces, dashed lines to the 
symmetric slab for the Ga-t surface. Macroscopic averages 
vanish in the bulk and in the vacuum, and their integrals 
over the surface regions give the surface stresses. Solid and 
dashed lines are distinguishable only in the region of the As-t 
surface. 



As already mentioned, for GaAs (111) there is no way 
to set up a symmetric slab containing only one kind of 
surface, since (111) and (111) surfaces are intrinsically 
different, and not related by any symmetry operation. 
We then we consider a comparison between the structure 
with both Ga-t and As-t (111) surfaces, and that one 
with Ga-t (111) and (ITT) surfaces (Fig. |). The solid 
lines represent planar and macroscopic averages of the 
slab with chemically different surfaces, the dashed ones 
are the analogous quantities for the Ga-t surface slab. 
Both the structures are now affected by an equal spuri- 
ous charge transfer ~ 0.04 electrons per cell; this is not 
unexpected since electron counting assigns in both cases 
the same occupation numbers to partially-filled surface 
bonds, that is, 3/4 on Ga-t (111) and 1/4 on As-t (111) 
and Ga-t (111). The comparison of the half-slab stresses 
containing the same Ga-t surface, i.e. the Ga-t surface 
stress, show an almost perfect agreement: we obtain - 
0.515 eV for the same-ion surface slab, and -0.508 eV 
for the different-ion surface slab. This is due both to 
the large number of bulk layers ensuring that there is no 
interaction between the surfaces in the slab, and to the 
(partially fortuitous) equality of charge transfer in the 
two slabs. 
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FIG. 3. Planar and macroscopic averages of stress density 
for GaAs (111). Solid lines refer to the slab with both Ga-t 
and As-t surfaces, dashed lines to that with (111) and (111) 
Ga-t surfaces. 



C. Ge/GaAs (100) and (110) interfaces 

In this Section we present calculations on Ge/GaAs 
interfaces. For interfaces, the computing effort saved by 
using space-resolved quantities can be even hig her th an 
for the surfaces. First, the considerations of Sec. [II B for 



polar surfaces apply in full also to inequivalcnt interfaces. 
Second, the simulation cell sizes are generally larger than 
for surfaces, hence the ability of avoiding separate calcu- 
lations is very much welcome. Third, the relevant quan- 
tities entangled in the integrated stress (or energy) are 
at least three, namely the interface excess stress and the 
bulk stress of the two bulk phases, so that more separate 
bulk calculations are needed within the usual approach 
than in the surface case. 

As an example of stress density application we first 
present the case of the non-polar abrupt lxl (110) 
Ge/GaAs interface. In Fig. [|, planar and macroscopic 
averages of the stress density are shown for a (GaAs)s 
(Ge)5 slab. The in- plane lattice constant is fixed at the 
theoretical value for bulk GaAs; since the equilibrium lat- 
tice constant of bulk Ge is slighty larger, the Ge film is 
under tensile stress. The structure is unrelaxed, and the 
interlayer distances are those of bulk GaAs. [Note that 
these interfaces are stoichiometric, hence Eq. ( |38|) can 
be used to get rid of the spurious bulk stress. The highly 
converged plane-wave basis used for the polar surfaces 
is not mandatory any more, and the cutoff was lowered 
to 25 Ryd. This produces an appreciable residual stress 
in GaAs.] The stress induced by the interface is barely 
visible in Fig. ^; only a small dip in the macroscopic 
average can be distinguished between the interface GaAs 
and Ge layers. For a non polar junction this is quite rea- 
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sonable, since the charge transfer across the interface is 
too small and interface-localized to produce a significant 
excess stress. Integrating the planar (or macroscopic) 
average over the bulk regions we obtain T^ s =-0.94 
eV/atom and T^f k =-1.52 eV/atom. The resulting in- 
terface stress of T mt = 0.02 eV/cell - 1 meV eV/A is 
practically vanishing. Clearly this nearly-matched inter- 
face is a test case where the stress due to the interface 
formation does not play any major role. Nevertheless, 
even in this case the power and semplicity of the stress 
density approach is evident. 
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FIG. 4. Planar average (solid line) and macroscopic aver- 
age (dashed line) of the stress density for the (110) Ge/GaAs 
interface. The presence of the interface is barely visible in the 
macroscopic average, since the interface is non polar and only 
very slightly mismatched. Only half slab is shown. 
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FIG. 5. Planar (solid line) and macroscopic (dashed line) 
average of the stress density for the mixed 2x1 (100) 
Ge/GaAs interface. GaAs and Ge are interfaced by a mixed 
Ge/Ga layer. By symmetry, only one interface is shown. 

The theoretical bulk GaAs in-plane lattice constant is 
imposed to the structure; the value of the macroscopic 
average in the bulk GaAs region represents just the finite- 
size basis error, whereas the bulk Ge region is also subject 
to strain. The bulk stresses are r^ As =-0.95 eV/atom, 
and T^ k =-1.61 eV (the small differences from the val- 
ues obtained in the (110) Ge/GaAs interface is due to 
k-point sampling). The interface region extends (stress- 
wise) over three layers, the mixed interface layer and 
the two adjacent to that. Outside this restricted region, 
the macroscopic average quickly recovers the equilibrium 
bulk value. The interface stress is 0.21 eV/atom, i.e. a 
condition of small compressive stress induced by the Ge 
epilayer. 



A more interesting case of interface-induced stress is 
illustrated in Fig. where we show planar and macro- 
scopic averages of the stress density for the mixed 2x1 
(100) Ge/GaAs interface. This polar interface is simu- 
lated with two bulk regions comprising 5 Ge and 5 GaAs 
layers, matched by a mixed Ga-Ge layer (12 layers and 
24 atoms in total). The slab contains an equal number 
of Ga and As atoms, so that the bulk GaAs contribution 
can be unambiguously subtracted out__ It is also easily 
verified by the electron-counting ruleE3 that no metal- 
lic states are present at the interface, which is in fact a 
reasonable candidate as the actual interface structure. 



D. Microscopic analysis based on the stress density 

So far we have been giving examples of applications of 
the stress density to the calculation of macroscopic stress 
in cases where direct calculations would be troublesome. 
Actually, however, the stress density appears to be most 
useful in analyzing processes driven or simply influenced 
by the stress, as a kind of "stress microscope". Recon- 
struction, adsorption, and growth are some of the phe- 
nomena in which stress can play a relevant role, possibly 
be a driving force, or simply function as a useful indi- 
cator of the ongoing processes. BI13 The integrated stress, 
though, is sometimes of limited use, and occasionally mis- 
leading. While it does tell us whether the application of 
strain will be energetically favorable, it easily misses the 
microscopic details, and cannot answer several legitimate 
questions such as: How far does the perturbation propa- 
gate below a surfacfi-pr across an interfaces? How is the 
stress redistributed^ into the perturbed region? -What 
interplay of tensile and compressive contributionsca pro- 
duces the resulting integrated stress ? Sometimes an- 
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swering these questions is fundamental to attain a con- 
sistent qualitative view of a specific process, and indeed 
the stress density has been usefully applied-ijO-this fash- 
ion to the study of surface reconstructions EjO Here, in 
particular, we discuss briefly the hex reconstruction of Ir 
(100), an apt example of what we mean by a microscopic 
analysis performed using stress density. 

In this reconstruction, a quasi-hexagonal buckled layer 
(whereby 6 atoms take place into a 5x1 surface area) 
is formed on top of the square (100) substrate. The in- 
crease in surface atomic density is expected to have a 
close relation to the surface stress, which is tensile and 
very large (-1.86 eV/lxl area) on the unreconstructed 
surface. The easy guess (the "stress hypothesis") would 
be that the densifying reconstruction is driven by the en- 
ergy gained by partially or completely relieving the stress 
of the pristine surface. This hypothesis was indeed made 
plausible by comparing the stress-driven energy gain with 
model estimates of the substrate-surface rebonding cost, 
based e.g. on bond cutting arguments, or mode ls of the 
defects needed to create the reconstruction.liStS 
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FIG. 6. Planar average of the stress density for unrecon- 
structed (dashed line) and 5x 1-reconstructed (solid line) Ir 
(100). The integrated stress is negative (i.e. tensile) for both 
the surfaces, and larger in magnitude for the reconstructed 



The validity of previous modelsll3 does not depend 
on the value of the stress after the reconstruction; it is 
reasonable to assume, though, that the stress should de- 
crease upon reconstructioiij-iittingly to its role of driving 
force. However, we founcOL-il that the surface tensile 
stress actually increases in magnitude by about 20% or 
0.3 eV/lxl area (on average on the x and y stress com- 
ponents of the anisotropic reconstructed surface) xm Ir 
(001). Similar results are suggested by experimentsE3 for 
Au (001). These new evidences on the integrated stress 
have increased the confusion, leading some to concludeO 
that the "stress hypothesis" is invalid, i.e. that stress 
is not the driving force to reconstruction, and enhancing 
the mismatch between data and physical intuition. It 
turns out that the application of stress densityEH provides 
a natural explanation of the reconstruction mechanism, 
reconciling all previous theoretical and experimental ob- 



servations. 

In Fig. ^| the stress density planar average is shown 
for the (fully relaxed) reconstructed and unreconstructed 
(100) surfaces. Upon reconstruction the outermost nega- 
tive and two outermost positive peaks, storing tensile and 
compressive contribution, respectively, are strongly re- 
duced. The surface layer relaxes outwards by a huge 20% 
of the ideal interplanar distance, which is also reflected 
in the outward shift of the tensile peak. The smearing of 
the latter is a token of the large surface buckling. From 
these findings, the following picture readily emerges. The 
decreased interatomic distance in the surface layer corre- 
lates with the severe reduction of the tensile stress within 
that layer (suppression of tensile peak). At the same 
time, the surface layer needs to relax outward to reduce 
its interaction with the substrate potential; this displace- 
ment correlates with the decrease of the compressive (ki- 
netic) contribution to the stress, because more room is 
provided to the electronic charge to delocalize in between 
the top surface layers. Thus the net negative change of 
total surface stress (more tensile after reconstruction) re- 
sults from the competing reductions of tensile and com- 
pressive contributions that are easily related to struc- 
tural changes. In particular, while the pre-existing ten- 
sile stress can be clearly identified as driving force, the 
change in compressive stress resulting from the surface 
rearrangement to accomodate the surface-substrate mis- 
match can be interpreted as a residual "resisting" force 
of the substrate. The post-reconstruction stress contains 
both contributions, inextricably entangled. It appear 
that the conclusion drawn from the post-reconstruction 
total stress, that the transition is not-stress driven, is 
incorrect. Further details on this reconstruction can be 
found elsewhere. lill 

The objection to this usage of the stress density based 
on its formal gauge dependence was considered in Sec. 
[II C| ; there we showed that perturbation-induced changes 
of different stress densities subject to the constraint 
of producing the same integrated stress are identical 
to linear order in the density, hence effectively gauge- 
independent. Thus, the analysis of stress density changes 
upon reconstruction (or adsorption, etc.) is indeed mean- 
ingful. 

It is worth remarking again here that the analysis of 
stress density changes upon reconstruction (or adsorp- 
tion, etc.) is meaningful, given the proof in Sec. [II C\ 
that perturbation-induced changes of different stress den- 
sities subject to the constraint of producing the same in- 
tegrated stress are identical to linear order in the density. 



IV. SUMMARY 



Drawing on the theory of quantum mechanical stress, 
we introduced the stress density in density functional the- 
ory, giving specific prescriptions for a practical and effi- 
cient implementation for the plane wave ultrasoft pseu- 
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dopotcntial method within the local-density approxima- 
tion. In analogy with the Chetty-Martin energy den- 
sity, the stress density provides a spatial resolution of 
the contributions to the integrated macroscopic stress 
tensor. We applied the newly introduced concept to a 
wide-ranging set of test cases selected from surface and 
interface physics. It is appropriate to point out that in 
this paper we gave just a flavor of the possible applica- 
tions, that are far from being exhausted by those pre- 
sented here. In perspective, the stress density is way 
more powerful and predictive than the integrated macro- 
scopic stress, and opens up new ways towards an under- 
standing of the physics of stress-driven phenomena. 
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TABLE I. Layer-by-layer (LBL) decomposition of the in- 
tegrated stress for an Al (111) surface slab. The slab is bulit 
up by 15 (9 atomic plus 6 vacuum) layers, labeled from -7 to 
+7. Layer is the bulk-like slab center, and +n and -n layers 
are identical by symmetry. Values are in eV/atom. 






±1 


±2 


±3 


±4 


±5 


±6 ±7 


r LBL -0.779 


-0.734 


-0.706 


-0.771 


-2.059 


0.729 


0.0 
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